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ABSTRACT 

A  linear  triangular  finite  element  formulation  was  used 
to  solve  the  steady  state  two  dimensional  conduction  heat 
transfer  equation.   A  FORTRAN  IV  computer  program  of  the 
above  formulation,  employing  double  precision  arithmetic 
and  compact  storage  techniques,  was  applied  to  study  heat 
transfer  in  an  internally  finned  rotating  heat  pipe.   Results 
were  obtained  for  water  in  a  copper  and  stainless  steel 
condenser  with  varying  outside  heat  transfer  coefficient, 
rotational  speed,  and  fin  angle  (number  of  fins) . 

Numerical  results  of  the  heat  transfer  rate  in  the  copper 
condenser  were  shown  to  have  less  than  2%  variance  to  those 
obtained  earlier.   A  rotating  heat  pipe  designed  to  operate 
with  a  small  value  of  outside  heat  transfer  coefficient 
experiences  no  significant  increase  of  the  heat  transfer  rate 
with  an  increase  in  rotational  speed,  since  the  value  of  the 
outside  thermal  resistance  dominates. 

Heat  transfer  rate  continuously  increases  as  the  fin 
angle  decreases.   However,  the  increase  is  only  slight  when 
the  fin  angle  is  less  than  11  degrees. 
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I.   INTRODUCTION 

A.   THE  ROTATING  HEAT  PIPE 

The  rotating  heat  pipe  is  a  closed  container  designed 
to  transfer  a  large  amount  of  heat  in  rotating  machinery. 
Its  three  main  component  parts  are:   a  cylindrical  evaporator, 
a  truncated  cone  condenser,  and  a  working  fluid  as  shown 
in  Figure  1. 

At  rotation  above  the  critical  speed  of  a  rotating  heat 
pipe,  the  working  fluid  forms  an  annulus  in  the  evaporator, 
and  will  be  vaporized  by  heat  addition  in  it.   The  vapor 
flows  toward  the  condenser,  as  a  result  of  a  pressure  differ- 
ence, transporting  the  latent  heat  of  evaporation  with  it. 
External  cooling  of  the  condenser  causes  the  vapor  on  the 
inner  wall  to  condense  and  release  its  latent  heat  of  evapor- 
ation.  The  centrifugal  force  due  to  the  rotation  has  a 
component  acting  along  the  condenser  wall  that  will  act  to 
drive  the  condensate  back  to  the  evaporator  where  the  cycle 
is  repeated. 

In  a  conventional  heat  pipe,  the  force  driving  the  con- 
densate back  to  the  evaporator  is  due  to  capillary  action, 
which  poses  a  limit  to  its  operation.   However,  the  driving 
force  in  the  rotating  heat  pipe  can  be  designed  to  operate 
in  any  orientation,  including  in  a  zero  gravity  field. 
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B.   UTILIZATION  OF  AN  INTERNALLY  FINNED  CONDENSER 

The  first  theoretical  investigation  of  the  rotating 
heat  pipe  conducted  at  the  Naval  Postgraduate  School,  was 
performed  by  Ballback  [1]  in  1969.   He  studied  the  limi- 
tations of  performance  imposed  on  the  rotating  heat  pipe 
due  to  various  fluid  dynamic  mechanisms.   Using  existing 
theory  and  experimental  correlations,  he  was  able  to  esti- 
mate the  sonic  limit,  boiling  limit,  entrainment  limit, 
and  the  condensing  limit  of  performance. 
1 .   The  Sonic  Limit 

When  increasing  the  heat  flux  in  a  rotating  heat 
pipe,  it  is  entirely  conceivable  to  reach  a  limiting  flow 
rate  of  the  vapor  brought  on  by  a  choked  flow  condition  in 
the  pipe.   This  condition  imposes  a  limiting  value  on  the 
amount  of  energy  the  vapor  can  transport,  thus  reducing  the 
effectiveness  of  the  heat  pipe.   Since  the  heat  transfer 
method  in  the  rotating  heat  pipe  depends  strongly  on  the 
latent  energy  of  the  working  fluid,  the  rate  of  heat  trans- 
fer in  the  heat  pipe  can  be  expressed  as: 


Qt  =   mv  h  CI-D 


where   m  =  vapor  mass  flow  rate  in  lbm/hr .   From  the  con- 


tinuity equation, 


mv   =   pv  Uv  A  (1-2 
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where 


U   =   velocity  of  the  vapor  in  ft/sec,  and 

2 
A   =   cross  sectional  area  for  the  vapor  flow  in  ft 


The  heat  transfer  rate  becomes 


Q4.   =   P   U   A  h.  (1-3) 

t       v   v     fg 


and  the  vapor  velocity  is  considered  to  be  sonic, 


Uv   =   c   =   /gQkRT  (1-4) 


where 


c    =   sonic  velocity  in  ft/sec 

2 

g0   =   gravitational  constant  =  32.1739  ft-lbm/lbf-sec 

k    =   ratio  of  specific  heats 

R   =   gas  constant  in  ft-lbf/lbm  R,  and 

T    =   absolute  temperature  in  °R. 

2 .   Boiling  Limit 

It  has  been  postulated  by  Kutateladze  [2]  that  the 
transition  from  nucleate  to  film  boiling  is  totally  a 
hydrodynamic  process.   With  this  postulation,  he  determined 
the  theoretical  formula  for  predicting  the  burnout  heat 
flux: 
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1/4 

Qt      =      K/pv"  Ab   hfg{a    g(pf    "    pv)}  (I"5) 


2 


where 


K    =   constant  value 

3 

p     =   density  of  the  vapor  in  lbm/ft 

A,    =   heat  transfer  area  in  the  boiler  in  ft 
hf    =   latent  heat  of  vaporization  in  Btu/lbm 

a    =   surface  tension  in  lbf/ft 

2 
g     =   acceleration  of  gravity  in  ft/hr 

pf    =   density  of  fluid  in  lbm/ft 

The  experimental  data  obtained  by  Kutateladze  suggested  a 
value  for  K  in  the  range  of  0.13  to  0.19. 
3 .   Entrainment  Limit 

The  flooding  constraint  in  a  wickless  heat  pipe  was 
examined  by  Sakhuja  [3],  who  assumed  that: 

a)  The  vapor  and  liquid  are  exposed  to  each  other 
in  a  counter flow  movement. 

b)  The  counterf lowing  vapor  tends  to  retard  the 
falling  film  of  liquid  through  shear  stress. 

c)  The  liquid  film  remains  stable  and  smooth,  and 
the  shear  stress  is  usually  small. 

d)  The  buoyancy  forces  resulting  from  the  density 
difference  between  the  vapor  and  liquid  are  the 
means  for  maintaining  the  counterf low. 
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e)   The  buoyancy  forces  are  balanced  by  dissipative 
effects  which  are  proportional  to  the  momentum 
fluxes  of  vapor  and  liquid  streams. 
His  resulting  correlation  for  flooding  is 


,2 


A  C   h,.  /gD(p--p  )p 
q    _   _x fg  ^   Mf  Hv  yv  (I_6 

fc       [1  +  (pv/pf)1/4]2 


where 


Q    =   heat  transfer  rate  in  Btu/hr 

2 
A    =   flow  area  in  ft 
x 

C    =   dimensionless  constant,  0.725  for  tube  with 
sharp  edged  flange 

hf    =   latent  heat  of  vaporization  in  Btu/lbm 

2 

g     =   acceleration  due  to  gravity  in  ft/hr 

D     =   inside  diameter  of  heat  pipe  in  ft 

3 

p.p   =   density  of  the  fluid  in  lbm/ft 

p     =   density  of  the  vapor  in  lbm/ft'"' 

4 .   Condensing  Limit 

Ballback  [1]  determined  the  condensation  solution 
for  a  rotating  heat  pipe  by  modeling  the  condenser  section 
of  a  rotating  heat  pipe  as  a  rotating  truncated  cone,  from 
which  the  following  condensation  limit  was  developed: 


~  k,.  p-  co   h^  [T  -T  ]   4  T      T  X 

rs  r2  ff      fgs   wir,_  _  .    ,,3    _  3  -i  4 

Q.   =   tt{t  ^-^ }  {  [R   +  Lsm  ($>  ]       -  R   } 

u  j  Z  O  O 


U.p  sin 


(1-7) 
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where 


\ 


Q.  =  total  heat  transfer  rate  in  Btu/hr 

kf  =  thermal  conductivity  of  the  condensate  film 
in  Btu/hr-ft-F 

3 

pf  =  density  of  fluid  in  lbm/ft 

co  =  angular  velocity  in  1/hr 

h-  =  latent  heat  of  vaporization  in  Btu/lbm 

T  =  saturation  temperature  in  F 

T  =  inside  wall  temperature  in  F 

yf  =  viscosity  of  fluid  in  lbm/ft-hr 

(J)  =  half  cone  angle  in  degrees 

R  =  minimum  wall  radius  in  ft 
o 

L  =  length  along  the  wall  of  the  condenser  in  ft 


The  condensing  limit  equation  is  a  function  of  the  geometry 
and  speed  of  the  rotating  heat  pipe,  and  the  physical 
properties  of  the  working  fluid. 

Tantrakul  [4]  calculated  these  limitations  on  a 
heat  pipe  with  specific  physical  characteristics  as  shown 
in  Table  1,  with  the  results  shown  in  Figure  2. 

TABLE  1 

Specification  of  a  Typical  Rotating  Heat  Pipe 

Length  14.000  inches 

Minimum  diameter      2.000  inches 

Wall  thickness      .  0.125  inches 

Internal  half  angle   1.000  degree 

Rotating  speed        2700  RPM 
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Obviously  from  the  results  in  Figure  2,  the  condensing  limit 
line  is  the  predominant  limitation  for  the  amount  of  heat 
that  can  be  transferred  from  the  heat  pipe.   However,  the 
other  limitations  may  become  important  as  the  heat  pipe 
geometry  and  operating  conditions  are  varied. 

In  order  to  augment  the  heat  transfer  capacity  of 
the  heat  pipe,  the  current  efforts  are  aimed  at  raising  the 
condensing  limit  line  which  may  be  accomplished  by: 

a)  a  high  value  of  cone  angle,  to  increase  the 
driving  force. 

b)  some  type  of  promoter  of  dropwise  condensation 
to  increase  the  value  of  the  convective  heat 
transfer  coefficient,  h. 

c)  use  of  an  internally  finned  condenser  to  increase 
the  inner  wall  surface  area  and  the  value  of 

h,  since  the  presence  of  the  fin  will  decrease 
the  condensate  film  thickness. 
By  machining  fins  in  the  inner  wall  of  the  condenser, 
a  significant  improvement  of  heat  transfer  will  be  realized. 

C.   ANALYSIS  OF  THE  INTERNALLY  FINNED  ROTATING  HEAT  PIPE 

A  disadvantage  of  improving  the  performance  of  a  rotating 
heat  pipe  with  fins  machined  in  the  inner  wall  is  that  an 
increase  in  the  wall  thermal  resistance  results  due  to  the 
increase  in  wall  thickness  brought  on  by  the  fins.   By  using 
a  material  with  a  high  value  of  thermal  conductivity,  this 
effect  can  be  minimized. 
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Figure  3   Internally  Finned  Condenser 

Geometry  Showing  Fins,  Troughs 
and  Lines  of  Symmetry 


22 


o^_ 


h 


out 


Figure  4   Internally  Finned  Condenser  Geometry 
Considered  in  Analysis  of  Schafer  [5] 
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Schafer  [5]  developed  an  analytical  model  for  the  case 
of  fins  with  a  triangular  profile  as  shown  in  Figures  3  and 
4.   In  the  analysis  he  assumed  one  dimensional  heat  conduc- 
tion through  the  wall.   Corley  [6]  for  this  same  case 
developed  a  two  dimensional  heat  conduction  model  using  a 
Finite  Element  Method,  and  also  assumed  a  parabolic  tem- 
perature distribution  along  the  fin  surface.   His  results 
indicated  a  significant  improvement  in  heat  transfer  per- 
formance of  about  75%,  above  that  predicted  by  the  one 
dimensional  model  of  Schafer  [5].   However,  Corley  [6] 
cautiously  noted  a  probable  error  of  50%  existed  in  the 
calculation  of  the  heat  transfer  at  the  fin- apex,  and  conse- 
quently mentioned  that  there  may  be  a  total  heat  transfer 
error  of  as  high  as  15%.   Also,  he  noted  that  his  analysis 
was  not  valid  fro  a  stainless  steel  condenser  because  of 
numerical  difficulties.   Tantrakul  [4]  modified  Corley' s 
computer  program  by  increasing  the  number  of  finite  elements 
in  order  to  minimize  the  apex  heat  transfer  error.   His 
results  with  this  modification  converged  with  the  results 
of  Corley. 

D.  THESIS  OBJECTIVES 

The  objectives  of  this  thesis  are: 

1.  To  generate  a  general  computer  program  for  two 
dimensional  steady  state  Conduction  Heat  Transfer. 

2.  To  use  this  program  to  study  heat  transfer  in  a 
finned  rotating  heat  pipe. 
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3.   To  conduct  a  parametric  study  for  a  proposed  heat 
pipe  design. 
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II.   FINITE  ELEMENT  SOLUTION 

A.   REVIEW  OF  THE  PREVIOUS  ANALYSIS 

Schafer  [5]  studied  the  one  dimensional  model  heat 
transfer  solution  and  Corley  [6]  studied  the  two  dimensional 
model  for  an  internally  finned  rotating  heat  pipe.   Both 
used  the  same  assumptions  and  boundary  conditions  based 
upon  the  analysis  of  Ballback  [1] ,  which  are  similar  to 
those  used  in  the  Nusselt  analysis  of  film  condensation  on 
a  flat  wall.   The  more  important  of  those  assumptions  are: 

1.  steady  state  operation 

2.  film  condensation,  as  opposed  to  dropwise  condensation 

3.  laminar  condensate  flow  along  both  the  fin  and  the 
trough 

4.  static  balance  of  forces  within  the  condensate 

5.  one  dimensional  conduction  heat  transfer  through 
the  film  thickness  (no  convective  heat  transfer  in 
the  condensate  film) 

6.  no  liquid  -  vapor  interfacial  shear  forces 

7.  no  condensate  subcooling 

8.  zero  heat  flux  boundary  conditions  on  both  sides  of 
the  wall  section  (symmetry  conditions),  as  shown 

in  Figures  5  and  6 

9.  saturation  temperature  at  the  fin  apex 

10.  zero  film  thickness  at  the  fin  apex 

11.  negligible  curvature  of  the  wall 
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The  complete  development  of  Schafer's  finned  condenser 
heat  transfer  is  presented  in  his  thesis.   His  results  are 
used  as  a  comparison  in  this  thesis. 

However,  since  the  same  fluid  flow  and  heat  transfer 
theory  applies  to  the  two  dimensional  conduction  model, 
Corley  [6]  utilized  this  model  and  further  developed  it 
using  the  Finite  Element  Method.   He  considered  an  internally 
finned  geometry  as  shown  in  Figure  6 ,  and  used  the  fin 
condensate  film  equation  derived  by  Schafer  [5] 


.  k  (T  -T  (z) )  y-  dz 

5  (z) d5  (z) 


2   2 

p-   co  (R,  +  x  sin  |)  h-   cos  $  cos  a 


(n-i: 


where 


6     =  the  fin  condensate  film  thickness  in  ft" 

kf   =  thermal  conductivity  of  working  fluid  in 
Btu/hr-ft-F 

3 

pf    =  density  of  working  fluid  in  lbm/ft 

yf    =  viscosity  of  working  fluid  in  lbm/ft-hr 

co     =  angular  velocity  in  1/hr 

Rn    =  minimum  radius  in  ft 

x    =  distance  along  apex  of  fin  in  ft 

Z     =  distance  along  the  fin,  from  the  apex  in  ft 

0     =  half  cone  angle  in  degrees 

hf    =  latent  heat  of  vaporization  in  3tu/lbm 

He  assumed  a  parabolic  temperature  distribution  along  the 
fin  surface, 
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Tw(z)   =   rxz2  +  Slz  +  t±  (II-2) 


to  calculate  the  convective__heat  transfer  coefficient  along 
the  fin  surface.   He  developed  a  solution  of  the  steady 
state  two  dimensional  conduction  heat  transfer  equation: 

Vk  •  VT   =   0  (II-3) 

with  boundary  conditions  as  shown  in  Figure  5 
a)   saturation  temperature  at  the  apex, 


T   =   T 

sat 


b)   convective  boundary  condition  in  region  [1] 


9T 
■k  ^—     =   h, (T  -  T   . ) ,   and  in  region  [3] 
dn       1       sat 


•k  |£  =   h0(T  -  t  ) . 
8n       2       °° 


3T 
c)   symmetry  conditions,  ^—  =  0 ,  in  region  [2]  and 

dn 

region  [4 ] . 

He  used  the  8  node  quadratic  isoparametric  element  in  the 

mesh  shown  in  Figure  6 ,  and  a  modification  of  the  three 

dimensional  isoparametric  Finite  Element  Method  computer 

program  assembled  by  Lew  [7] . 

Applying  the  boundary  condition  a)  ,  T  =  T   ,  at  z  =  0 

sat- 


gives  t,  =  T     and  equation  ( 
1     sat      ^ 


II-2)  becomes: 
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2 

Tw(z)   =   rxz   +  sl2  +  Tsat,  (II-4) 


or 


Tsat  "  Tw(2)   "   "rlz2  "  slz'  (I1-5) 

Upon  substituting  equation  (II-5)  into  equation  (II-21, 
and  integrating  the  result  from  z  =  0  to  z ,  gives  an 

equation  for  6(z): 

■ 

z3        2 
4  k  (-ri  \   -  s   V)uf 1/4 

5(z)   =   [ — ^ 2 ' 

pf   co  (Rfl  +  x  sin  <£)  hf   cos  $  cos  a 


»  • 


(II-6) 


Assuming  one  dimensional  conduction  heat  transfer  through 
the  thin  condensate  film,  the  local  convective  heat  transfer 
coefficient  is: 


3         2       2 

k_  k-      p_      co    (Rrt    +  xsin    p)h£    cos   cbcos   a    1/4 

h(z)      =   —     =       r    f    f - ^2 ] 

niz;fin      o(z)      L  z3       z2 

4  yf(-rx  —  -  s±   -y) 

(II-7) 


where   kf  =  fluid  thermal  conductivity  in  Btu/hr-ft-F  which 
is  applicable  from  z  =  0  to  z  =  Z*  =  Z»  -  6*/cos  a  ;  where 
6*  is  the  trough  condensate  film  thickness,  as  shown  in 
Figure  5.   Across  the  width  of  the  trough,  the  convective 
heat  transfer  coefficient  is: 
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kf 

h(x)  trough   =  WWf  (II-8) 


The  heat  transfer  rate  for  the  entire  condenser  is  then 
determined  within  the  computer  program  algorithm  by  dividing 
its  axial  length  into  100  increments  of  length  Ax  (see 
Figure  4),  and  solving  for  the  incremental  heat  transfer 
rate  using  an  iterative  procedure.   The  incremental  rates 
are  then  summed  to  yield  the  total  heat  transfer  rate. 

Figure  7  shows  a  flow  chart  description  of  the  computer 

program  used  in  Corley's  analysis  [6].   Constant  values  of 

r,  and  s,  are  determined  at  each  increment  using  Lagrangian 

interpolation  fit  of  the  temperature  at  the  nodes  1,  3, 

and  6  (see  Figure  6)  from  the  previous  iteration.   At  the 

first  increment,  values  of  r,  and  s,  are  calculated  from  an 

initial  guess  of  the  temperature  at  these  nodes.   The  average 

nodal  values  of  h,    ,  h_.     and  the  average  of  h  of  the 

lavg    3avg 

left  hand  boundary  region  associated  with  node  6  (see  Figure 
6),  are  calculated  using  the  position  dependent  value  h(z)--;n 
(equation  (II-7) )  in  a  Simpson's  Rule  approximation.   In 
the  trough,  h(x)      .  is  calculated  from  equation  (II-3), 
using  the  6*  value  calculated  from  the  previous  increment, 
and  this  value  is  used  as  the  average  convective  heat  trans- 
fer coefficient  for  nodes  9  and  12,  and  the  right  hand 
boundary  of  node  6.   The  average  heat  transfer  coefficients 
calculated  for  the  left-  and  right-hand  boundary  region  of 
node  6  are  weighted  by  their  respective  region  length.   These 
values  are  then  added  to  yield  the  average  convective 
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Dimensional  Conduction  Analysis  of 
Corley  [6] 
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coefficient  for  node  6 ,  hg    .   The  remaining  boundary 
conditions  are  now  applied,  and  the  incremental  heat  trans- 
fer rate  per  unit  of  condenser  length  Q.,  and  the  nodal 
temperature  values  are  generated  from  the  Finite  Element 
Method  subroutine.   If  the  resultant  heat  transfer  rate 
converges  to' within  a  0.05%  change  from  the  value  of  the 
previous  iteration,  i.e. 


Qj  -  Qi-i 

Q3 


<      0.0005  (II-9) 


then  the  increment  is  deemed  solved  and  Q.  is  set  equal  to 
the  final  iteration,  Q..   However,  if  the  convergence 
criterion  is  not  met,  newly  calculated  nodal  temperature 
values  are  used  to  establish  a  new  set  of  inner  wall  convec- 
tive  heat  transfer  coefficients,  and  the  Finite  Element 
solution  is  repeated  until  the  heat  transfer  rate  converges 
to  within  0.0005  tolerance.   The  incremental  mass  flow  rate 
in  the  trough  is  now  calculated  using  the  relation 

Q. 

M.     .   =  2  r^-Ax  (11-10) 

tot  i       hc 

fg 

where  this  value  is  added  to  the  mass  flow  from  the  previous 
increments.   A  modification  of  the  trough  mass  flow  rate 
equation  from  Schafer's  thesis  [5]  is  used  to  calculate  the 
subsequent  interval's  trough  condensate  thickness  5*(x), 
with  the  polynomial  root  finder  subroutine : 
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2   2  2 

p   oj  (Rn  +  xsin  <j))5*(x)   sin  <j>  „ 

Mtot(x)  = 3 [6*(x)e  +  6*(.x)  tan  a} 

(II-ll) 

This  incrementation  is  continued  until  the  entire  length  of 
the  condenser  is  transversed.   Incremental  heat  transfer 
rates  are  summed,  hence  the  total  heat  transfer  rate: 

100 

Q*.  4.   =   2  N-.    I  (Q.Ax)  (11-12) 

tot        fin  L         i 

i=l 

where  N _ .   is  the  number  of  fins  along  the  inner  wall 

fin  3 

circumference . 

To  start  the  incrementation,  two  initial  guesses  are 
made.   These  are  the  initial  value  of  5*  which  is  taken 
from  the  analysis  by  Sparrow  and  Gregg  for  condensation  on 
a  rotating  disk  [8]  and  the  initial  values  of  the  tempera- 
tures at  the  nodes  1,  3  and  6. 

Corley  [6]  noted  that  an  error  of  50%  due  to  calculation 
of  heat  transfer  rate  at  the  apex  may  cause  a  total  error 
of  15%.   To  minimize  this  error  he  suggested  an  increase  in 
the  number  of  finite  elements. 

Tantrakul  [4]  continued  Corley 's  study.  He  increased 
the  number  of  finite  elements  between  the  apex  of  the  fin 
and  the  trough  as  shown  in  Figures  8  and  9 .  To  save  com- 
puter execution  time,  he  reduced  the  number  of  increments 
to  25  instead  of  100  increments  as  in  Corley' s  calculation. 
However,  even  with  the  reduced  number  of  increments,  the 
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results  of  Tantrakul  [4]  converged  favorably  with  Corley's 
results  [6]  . 

B.   FORMULATION 

In  this  thesis  the  author  used  the  same  physical  rela- 
tion and  technique  to  calculate  the  heat  transfer  rate  from 
the  condenser.   These  are  listed  below: 

a)  assume  parabolic  temperature  distribution  along  the 
fin  boundary 

b)  utilize  equation  (II-6)  to  determine  the  condensate 
thickness ,  6  (z) 

c)  use  equation  (II-9) ,  (11-10) ,  and  (11-11)  to 
calculate  the  mass  flow  rate  in  the  trough,  and 
determine  the  subsequent  condensate  thickness. 

d)  iterate  to  arrive  at  the  desired  solution  with  the 
convergence  criterion  0.01,  since  the  results  are 
converged  at  this  value  (see  Table  2). 

For  the  Finite  Element  Method  solution,  the  author  used  a 
linear  triangular  finite  element  model  as  shown  in  Figures 
10  through  12. 

Corley  [6]  solved  the  field  equation  for  the  fin  with 
the  following  boundary  conditions : 

a)  temperature  at  the  apex  is  equal  to  the  saturation 
temperature 

b)  convective  boundary  conditions  in  the  inside  and 
outside  of  the  condenser 

c)  symmetry  condition  on  the  sides. 
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with  25  Linear  Triangular 
Finite  Elements 
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Figure  11   Condenser  Geometry  Considered 
with  4  0  Linear  Triangular 
Finite  Elements 
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Figure  12   Condenser  Geometry  Considered 
with  108  Linear  Triangular 
Finite  Elements 


42 


The  first  boundary  condition  indicates  no  heat  transfer 
at  the  apex  and  a  value  of  infinity  for  the  convective  heat 
transfer  coefficient  h.   Since  the  value  of  h  in  reality 
never  equals  infinity,  the  author  let  the  value  of  tempera- 
ture at  the  apex  float,  which  is  in  agreement  with  the 
theoretical  analysis  for  the  extended  surface  in  a  one 
dimensional  case.   Only  for  a  fin  with  a  small  angle  or  a 
very  long  length  can  the  temperature  at  the  apex  become 
equal  to  the  environment  temperature  (saturation  temperature 
in  this  case) . 

Therefore  the  statement  of  the  problem  for  the  formula- 
tion of  the  Finite  Element  Method  is  shown  in  Figure  13: 


2      2 

+  i-I  =   0  (11-13) 


i  2    ^  2 
3x     3y 


with  the  boundary  conditions 


a)  in  region  1,   -k  |£  =  h, (T  -' T      J         (11-14) 

3  an    1      sat 

b)  in  region  3,   -k  ^  =  hJT  -  T  )  (11-15) 

3  an    2      °° 

3T 

c)  in  region  2  and  4,    ^—  =  0  (11-16) 

^  3n 


Dividing  the  domain  into  linear  triangular  finite  elements , 
the  differential  equation  within  each  element  is: 


2  (e)    a2  (e) 

i-T  +  °  T      =   0  (11-17 

3x        3y 
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Figure  13   Differential  Equation  and  Boundary 
Conditions  Considered  in  Analysis 
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with  the  boundary  conditions: 


.,  ■-*(•>  |li!L   =   h'e»(T^-Ts(e))         (11-13 


for  elements  with  sides  along  the  convective  boundary,  and 
b)   §1—   =   0  (11-19) 

for  others.   In  the  above  equations, 

(e) 

T      =   temperature  within  the  element 

(e) 
k      =   thermal  conductivity  of  the  element;  assumed 

to  be  constant 

(e) 

h      =   convective  heat  transfer  coefficient  along 

the  boundary  of  the  element,  assumed  to 
be  constant 

(e) 

T      =   surrounding  temperature  at  the  boundary  of 

the  element. 

Define  the  approximate  value  of  the  temperature  distribution 
within  each  element  as : 

3 

T(e)   =   TN(e)  Cx,y)   =    I    Ni(x,y)  T.(e) 

1=1  (11-20) 

where 

N.      =   Two  dimensional  linear  shape  function,  and 

(e) 
T.      =   nodal  temperature. 

Using  equation  (11-20)  ,  equation  (11-17)  can  be  written  in 
approximate  form  as 
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3    fol  32N-     3    /  *  92N. 
R   =    [T.  (8)  — JL  +   Y  T.  (e)  — -i      (H-21 

i-1       9X      i-1  "     ^ 


where 


R      =   residual. 


Applying  the  Galerkin  criterion 


/R  N.  dA(e)   =   0  (11-22) 


and  substituting  equation  (11-21)  into  equation  (11-22), 
gives  : 


2  2 

d    N\        ,  ,  ,  ,  3  N.       ,  , 

)l   i    3*    J  (e)       ay   D 

e)  A[e) 


Applying  Gauss '  theorem  to  the  factors  on  the  right  side 
yields  the  symmetric  equation: 

,  ,        9N.9N.     3N.3N.     ,  . 

-7t  (e)      f(   x   J   +  _i__l)dA(e) 
*Ti        j(  3x  3x   +  3y  3y  } dA 

A(S) 

(e) 


ST.. 

+      (N.  — i dLie;  =   0      (11-23 

;  1    3n 


Lle 
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where 


(e) 
A     =   area  domain,  and 

(e) 
L     =   length  domain. 


Substituting  equation  (11-18)  and  (11-19)  into  the  second 
term  of  equation  (11-23),  gives: 


9T  (e) 

/N.  -4— -dL(e)   =   -      /N.  ^(T(S))  dL(e) 
'    j        an  J    j  _ 


(11-24) 


L(e)  L(e) 


where 


/  .  (T.  '    -  T  (  ')         for  an  element  with  its 
k  side  along  the  convec- 

tive  boundary 


0  others 


Therefore  equation  (11-23)  can  be  written  in  the  form 


3N.3N.    9N.3N. 


I    T.(e)      /[^-M-  +  t-V1]  dA(S) 

L    1  ;   dX  dx  dV  dV 


i-1        (e) 


+       [N.  ^>(T(e))  dL(e)   =   0       (11-25) 
L(e» 


47 


The  solution  of 


8N.SN.     aN,3N^ 
_ J j- 


(e) 


J-   ] 
3x  3x    '  3y  3y 


dA 


(e 


(pij)(S)    (11-26) 


where 


ID 


b.b. 
J-  3 

4A(6)    4Ate 


+ 


c  .  c  . 
i  3 


(11-27) 


b.    =   Y.  -  Y. 
i       j  k 


(11-28) 


c±        »   Xk  -  X.  (11-29) 

j  }   =   1,  2,  3 

k  ' 

as  shown  in  Figure  14.   Recalling  equation  (11-24),  for 
elements  with  an  edge  along  a  convective  boundary , 


!e) 


/N.0(T(e))dL(e)   =   T.(e)  & 
J    3  1  ^ 


(e) 


(N.N.  dL 


(e! 


!e) 


-  T 


(e)  h 


(e) 


TeT 


(e) 


/N.  dL 

D 


(e) 


(11-30) 


Applying  equation  (11-30)  for  elements  with  the  nodal  point 
1  and  2  located  on  the  convective  boundary,  gives 
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Figure  14   Linear  Triangular  Finite  Element  Geometry 
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(e)      yPx)<Pl'P2' 


X*T.<        -,     rp  f).   (  ,T  ( 


/N   ♦(Tw)dLw   =   T±w  2^      /  p2  dL 


L(S)  k     L(' 


(e)       \pl 


_  T  <e>  h,  r)0  ( dL(e) 

L(*>  (o) 

(11-31 


where 


p,,P2   =   one  dimensional  linear  shape  function 


Therefore,  equation  (11-23)  can  be  written 


<Pi  P9  o> 


lTi(e)[(pij)(e)  +^  j         P2( 

i  k     (e)   (  0 


Pi 

^S-T  (S)       /  ]P2(   dL(' 
Tie)   s         j  / 

L(e)   (  0 


or 


(K)(e){T}(e)   =   (F}(e)  (H-31) 
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1.1 


} 


4  3 


where 


A 

(K) 


(e)  _ 


/    . 

b  2+c  2  h(e)L   (e) 

4A(6)      3k(e) 

v242/e,^(e) 

4A(e)       6k(e) 

b  2^  2  .(e)    (e) 
b2  ^2  ,  h   L21 

/ 

b.b^+c.c- 
1  j  13 

y 

b2b3+C2c3 

4A(e)      6k(e) 
/ 

4A(e)     3k(e) 

/ 

b2b3+c2c3 

4A(^ 
b32+c32 

4A 


(e) 


4A 


(e) 


4A 


(e) 


and 


3 


{F} 


(e) 


h(e)  T  (e)     (e) 
s 21 

u(e) 


33    3)y  32 


The  element  matrix  equation  is  then  assembled  into  a  system 
matrix  equation  with  the  sequencing  based  upon  element  and 
system  nodal  point  connectivity. 

The  heat  transfer  rate  within  each  element  was  calculated 
by 


T  (e)    T  (e) 

Q(S)   -   h(e)  L12(S)  *x(-i 2 2- Ts(e) 


(11-32) 


where 


(e) 


(e) 


=   temperature  at  node  1 


T2  v     =   temperature  at  node  2 


9  O  9 
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(e) 
T      =   surrounding  temperature 

(e) 
h      =   convective  heat  transfer  coefficient 

(e) 
L»,     =   length  between  nodal  points  1  and  2 

as  shown  in  Figure  15. 

C.   THE  COMPUTER  PROGRAM 

The  program  consists  of  a  main  program  and  three 
subroutines ; 

a)  the  routine  "COORD"  used  to  define  positions  of 
system  coordinate  points 

b)  the  routine  "FORMAT "  used  to  formulate  the  Finite 
Element  Method  equations 

c)  the  routine  "BANDEC"  as  an  equation  solver  for  a 
symmetric  matrix  which  has  been  transformed  into 
banded  form. 

The  input  data  is  entered  with  the  units  of  inches  for  the 

2 

-length,  degrees  for  the  angle,  Btu/hr-ft  -F  for  convective 

heat  transfer  coefficients,  and  in  F  for  the  temperatures. 

Input  data  is  arranged  in  6  classifications,  and  is 
shown  in  Table  3  . 
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Figure  15  -  Heat  Transfer  Rate  of  an  Element,  Q 
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TABLE  3 
Classification  Of  Input  Data 

Classification  Type  of  input  data 

I  system  and  element  connectivity 

II  condenser  geometry 

III  statement  of  the  problem 

IV  convergence  criterion 

V  specification  to  determine  the  system 

coordinate  point 

VI  initial  guess  of  the  temperature  at  the 

nodal  points  along  the  fin  and  the  trough 

A  detailed  description  of  each  of  these  classifications 
follows : 

Classification  I 
READ  (5,  100)  NEL,  NSNP,  NBAN 
100  FORMAT  (315) 
Explanation:     NEL    =  number  of  elements 

NSNP   =  number  of  system  nodal  points 
NBAN   =  system  band  width 
READ  (5,  200)   (IEL, (ICOR (IEL, I) , 1=1 , 3) , IEL=1 ,NEL) 
200  FORMAT  (.415) 
Explanation:     IEL    =  the  element  number 

ICOR (IEL, I)    =  system  nodal  point  (NP)  corresponding 

to  NP  I  of  element  IEL 
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Note:   -  numbering  the  element,  first  for  all  elements  with 
a  convective  boundary,  starting  from  the  top  to 
the  bottom,  then  others 

-  numbering  the  nodal  point  of  an  element  is 
counterclockwise 

-  the  elements  with  convective-  boundary  have  nodal 
point   1   and   2   located  on  the  convective  boundary 

Classification  II 

READ  (5,  300)  CL,  CANGL,  RBASE,  R2 ,  THICK,  BFIN 
300  FORMAT  (  6G10.5) 
Explanation:      CL   =  condenser  length  in  inches 
CANGL   =  half  cone  angle  in  degrees 
RBASE   =  minimum  radius  of  the  condenser  in 
inches 

R2  =  RBASE-+  bkh 

THICK   =  condenser  wall  thickness  in  inches 
BFIN   =  the  height  of  the  fin  in  inches 
READ  (5,  40  0)  NDIV,  NEST,  NEFB ,  NBOTI ,  NBOTF 
400  FORMAT  (  515) 

Explanation:    NDIV   =  number  of  increments  along  the  length 

of  the  condenser 
NEST   =  number  of  the  element  on  the  right 

end  of  the  trough 
NEFB   =  the  element  number  with  convective 
boundary  located  at  the  base  of  the 
fin 
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NBOTI   =  the  element  number  with  convective 

boundary  located  at  the  right  hand 

of  the  bottom  side 

NBOTF   =  the  element  number  with  convective 

boundary  located  at  the  left  hand 

of  the  bottom  side 

Classification  III 

READ  (5,  500).  (MRPM(I) , 1=1 , 3) 

500  FORMAT  (  315  ) 

Explanation:     MRPM  =  number  of  RPM         (u>) 

READ  (5,  600)  (TSS ( I) , 1=1 , 3) 

600  FORMAT  (  3G10.5  ) 

Explanation:      TSS   =  saturation  temperature  of  the 

working  fluid         (T    ) 

READ  (5,  700)  TINF 

700  FORMAT  (G10.5) 

Explanation:     TINF   =  outside  temperature   (Tot) 

READ  (5,  800)   (HINF ( I ) , 1=1 , 3 ) 

800  FORMAT  (  3G10.5  ) 

Explanation     HINF   =  outside  convective  heat  transfer 

coefficient  (h   ,  ) 

out 

Classification  IV 

READ  (5,  9  00)  CRIT 
900  FORMAT  (  G10.5) 
Explanation:     CRIT   =  convergence  criterion  (0,01) 
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Classification  V 

READ  (5,  1000)   (FANGL(I) , (=1,3) 
1000  FORMAT  (  3G10. 5  ) 
Explanation:      FANGL  =  fin  half  angle       (a) 

READ  (5,  1100)  IFF 
1100  FORMAT  (  15  ) 
Explanation:       ,  IFF  =  n-1,  where  n  is  the  number  of  rows 

of  the  upper  triangular  fin  section 
READ  (5,  1200)   (KFIN (I) ,KFF (I ) , (=1 , IFF) 
1200  FORMAT  (  1615  ) 

Explanation       KFIN  =  the  number  of  system  nodal  points 

located  on  the  symmetric  boundary 
of  triangular  fin  section,  but  does 


not  include  the  system  nodal  points 
located  at  the  base  of  the  fin  and 
the  apex 
KFF  =  the  number  of  system  nodal  points 
located  along  the  fin  convective 
boundary,  but  does  not  include  the 
system  nodal  points  located  at 
the  base  of  the  fin  and  the  apex 
READ  (5,  1300)  DOBF ,  DOTH,  JTC,  JLC ,  JINT,  KT 
1300  FORMAT  (  2G10.5,  415  ) 
Explanation:       DOBF  =  number  of  column  within  the  fin 

DOTH  =  number  of  column  within  the  trough 
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JTC  =  number  of  the  system  nodal  point 
located  at  the  junction  of  the 
symmetry  boundary  and  the  line  of 
intersection  between  the  fin  and 
the  condenser  wall 
JLC  =  the  number  of  the  system  nodal 
point  located  at  the  center  of 
system  coordinate 
JINT  =  the  numerical  difference  between  the 
two  adjacent  system  nodal  points 
vertically  at  the  condenser  wall 
section 
KT  =  the  number  of  rows  within  the  wall 
section 
Note  :  as  an  example  for  25  finite  elements  is  shown  in 
Figure  16 
Classification  VI 
READ  (5,  1700)  T(IE) 
READ  (5,  1700)  T(IG) 
1700  FORMAT  (G10.5) 

Explanation:  T  =  initial  values  of  the  temperature 

estimate  at  the  nodal  points  1  and 
2,  for  elements  with  convective 
boundary  along  the  fin  and  the 
trough 
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IFF  =  2  DOBF  =3.0       JLC  =  17 

KFIN  =  2  ;  k    ^  DOTH  =  1.0       JINT  =  5 

KFF  =  3  •  6       JTC   =  7         KT  =  2 


Figure  16   Specification  for  Input  Data  to 
Determine  the  Coordinate  of  the 
System  Nodal  Points 
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D.   TEST  CASES 

Before  the  Finite  Element  Method  program  (routine 
"FORMAT")  was  applied  to  the  condenser  problem,  a  series  of 
test  cases  were  analysed.   These  cases  allowed  a  comparison 
to  the  results  obtained  from  previous  analyses  and  reference 
9.   These  cases  are  listed  below: 

1.  A  rectangular  block  in  a  uniformly  convective 
environment  on  three  sides,  with  a  specified 
temperature  on  the  non-convective  boundary.   These 
results  are  compared  with  reference  9  as  shown  in 
Figure  17. 

2.  A  rectangular  block  in  a  uniformly  convective 
environment  on  the  top  and  bottom  side,  insulated 
on  the  other  sides.   A  comparison  of  these  results 
with  the  exact  values  and  Corley's  results  [6]  is 
shown  in  Figure  18. 

3.  A  rectangular  block  with  triangular  fin  in  a  uniformly 
convective  environment  along  the  fin  and  bottom  side, 
insulated  on  the  wall.   A  comparison  of  the  results 
with  Corley's  results  is  shown  in  Figure  19. 

The  results  show  good  agreement  with  the  results  of  reference 
9  (case  1) ,  and  are  comparable  with  the  exact  values  (case  2) . 
Therefore,  "FORMAF"  can  be  used  for  further  applications. 
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-{ 


29894.  Btu/hr 
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Note:   numbers  in  parenthesis  refer  to  the  results  of 
reference  9 


Figure  17   Statement  and  Results  of  Test  Case  1 
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T    =    100°F 


h   =    100.    Btu/hr-ft    -deg   F 
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'  (7527.  Btu/hr) 


Note:   Numbers  in  parentheses  refer  to  results  of 
Corley  [6] . 


Figure  18.   Statement  and  Results  of  Test  Case  2 
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T,    =    200.     °F 


h1   =    100.    Btu/hr-ft2-deg   F  " 


0.75-C 


ft 


I 

V      0.5  ft 

1- 


T   =  100.  dF 
2 

h2  =  1000.  Btu/hr-ft2-deg  F 


c 

(ft) 

HEAT  TRANSFER  RATE  (Btu/hr) 

ONE      *) 
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DIMENSIONAL 

TWO 
DIMENSIONAL 

0.25 

6526. 

6727. 

6826. 

0.05 

6561. 

6754. 

6759. 

*)   Results  of  Corley  [6] 


Figure  19   Statement  and  Results  of  Test  Case  3 
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III.   NUMERICAL  RESULTS  AND  DISCUSSION 

The  analysis  concerns  itself  with  the  geometry  of  both 
a  copper  and  stainless  steel  condenser  as  tabulated  in 
Table  4,  using  water  as  the  working  -fluid. 

TABLE  4 

Specification  of  a  Typical  Internally 
Finned  Rotating  Heat  Pipe 

Condenser  length  =  9.000  inches 

Minimum  radius  =  0.775  inches 

Wall  thickness  =  0.03125  inches 

Height  of  fin  =  0.025  inches 

Cone  half  angle  =  1.0  degree 

Results  are  given  in  Tables  5  to  18  and  Figures  20  to  28. 

It  can  be  seen  from  the  data  that  the  solution  of  the  Finite 

Element  Method  converges.   This  can  be  seen  by  examining 

data  of  temperature  distribution  within  the  fins  of  25, 

40,  and  108  linear  triangular  finite  elements  (Tables  7  to 

9),  for  the  same  location  (nodes  10,  15,  and  28  respectively). 

The  heat  transfer  rate  solution  is  generally  affected 

in  the  same  manner  as  the  analysis  of  Schafer  [5]  when  the 

outside  convective  heat  transfer  coefficient  (h   ,),  rotational 

out 

speed  (oj)  ,  fin  angle  (2a)  ,  and  thermal  conductivity  (k   j,  ) 
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TABLE    9 


Theoretical    Heat   Transfer    Rate   on   a   Stainless    Steel   Finned 
Condenser    at    1000    RPM,    for   hQut   =    1000.       Btu/hr-ft   -deg.F 


Fin 

half 

angle 


10 


30 


50 


No. 

of 

fins 


276 


84 


40 


sat 
(°F) 


100 
150 
200 
100 
150 
200 
100 
150 
200 


25    EL 


8948.2 

23764.0 

38557.0 

8515.5 

22511.0 

36467.0 

8186.3 

21574.0 

34916.0 


Q(Btu/hr) 


40   EL 


8944.1 

23727.0 

38492.0 

8494.4 

22445.0 

36355.0 

8164.6 

21508.0 

34805.0 


108   EL 


8893.9 

23607.0 

38286.0 

8461.8 

22354.0 

36204.0 

8141.6 

21442.0 

34693.0 
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TABLE  10 


Theoretical  Heat  Transfer  Rate  on  a  Stainless  Steel  Finned 

Condenser  at  3000  RPM,  for  h     =  1000.  Btu/hr-f t2-deq.F 

out  3 


Fin 
half 

angle 

No. 

of 

fins 

T 
sat 

(°F) 

Q(Btu/hr) 

25  EL 

40  EL 

108  EL 

10 

276 

100. 

9096.4 

9100.4 

9058.1 

150. 

24193.0 

24200.0 

24098.0 

200. 

39277.0 

39284.0 

39123.0 

30 

■  84 

100: 

8818.4  ' 

8810.3 

8786.0 

150. 

23389.0 

23376.0 

23301.0 

200. 

37934.0 

37900.0 

37781.0 

50. 

40 

100. 

8584.1 

8580.7 

8561.1 

150. 

22720.0 

22690.0 

22642.0 

200. 

36823.0 

36768.0 

36689.0 
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TABLE  11 


Theoretical  Heat  Transfer  Rate  on  a  Stainless  Steel  Finned 
Condenser  at  1000  RPM,  for  h    =  5000.   Btu/hr-f t2-deg .F 


Fin 
half 

No. 

of 

fins 

T 
sat 

(°F) 

3(Btu/hr) 

angle 

25  EL 

40  EL 

108  EL 

10 

276 

100. 

26899.0 

26775.0 

26308.  o' 

150. 

70520.0 

68673.0 

200. 

113950.0 

113300.0 

110780.0 

30 

84 

100. 

22742.0 

22540.0 

22286.0 

150. 

58879.0 

58288.0 

57582.0 

' 

200. 

94734.0 

93748.0 

92581.0 

SO 

40 

100. 

20426.0 

20253.0 

20086.0 

150. 

52480.0. 

51978.0 

51513.0 

200. 

84275.0 

83443.0 

82676.0 
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TABLE  12 


i  u^a+-  Tr.nqfer  Rate  on  a  Stainless  Steel  Finned 
Theoretical  Heat  Transfer  Rate  o=        Btu/hr-f t2-deg .F 

Condenser  at  300U  KPM,  ror  uout 
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TABLE  13 

Temperature  Distribution  Within  the  Fin  at  the 

lot*1  Increment. 
(Copper  Condenser  with  25  Finite  Elements 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

98.99 

15 

98.08 

2 

98.79 

16 

97.94 

3 

98.77 

17 

98.00 

4 

98.60 

18 

97.99 

5 

98.58 

19 

97.97 

6 

98.53 

20 

97.93 

7 

98.43 

21 

97.80 

8 

9  8.41 

9 

98.35 

10 

98.23 

11 

98.01 

12 

98.18 

13 

98.17 

14 

98.14 

a  =  5 

3  deg 

T^  =  70°F 

h   4.  : 
on  t~ 

=  1000  Btu/hr-ft2 

-deg  F    RPM  =  10  00 

Tsat  : 

=  100°F 

Q  =  9697.2  Btu/hr 
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TABLE  14 

Temperature  Distribution  Within  the  Fin  at  the 

10tn  Increment 
(Copper  Condenser  with  40  Finite  Elements) 


NP 

T(°F) 

NP 

T  (  °  F ) 

NP 

T(°F) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

99.01 
98.84 
98.83 
98.69 
98.68 
98.64 
98.55 
98.54 
98.50 
98.45 
98.42 
98.41 
98.38 
98.32 

15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
2  8 

98.22 
98.04 
97.99 
98.17 
98.16 
98.14 
98.11 
98.06 
97.96 
97.92 
97.99 
97.98 
97.96 
97.94 

29 
30 
31 

97.90 
97.82 
97.79 

a  =  50  deg                      T   =  70°F 

J                                                           CO 

h     =  1000  Btu/hr-ft2-deg  F    RPM  =  1000 
out 

T     =  100°F                    Q  =  9687.2  Btu/hr 
sat 
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TABLE  15 

Temperature  Distribution  Within  the  Fin  at  the 

10th  increment 
(Copper  Condenser  with  108  Finite  Elements 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

99.04 

24 

98.39 

50 

97.90 

2 

98.91 

26 

98.33 

51 

98.05 

3 

98.90 

28 

98.20 

52 

97.96 

4 

98.79 

29 

98.05 

58 

97.85 

6 

98.76 

30 

97.99 

60 

97.84 

7 

98.68 

31 

97.96 

61 

97.97 

8 

98.68 

32 

98.27 

62 

97.96 

10 

98.63 

34 

98.25 

63 

97.95 

11 

98.58 

36 

98.20 

64 

97.94 

13 

98.56 

38 

98.12 

65 

97.92 

15 

98.51 

39 

98.02 

66 

97.90 

16 

98.49 

40 

97.96 

67 

97.90 

18 

98.47 

41 

97.94 

68 

97.88 

19 

98.45 

42 

98.14 

69 

97.82 

21 

98.37 

44 

98.04 

70 

97.78 

22 

98.41 

48 

97.91 

71 

97.76 

a  =  5 

0  deg 

Tc 

=  70°F 

h     4- 

out 

=  1000  Btu/hr-ft2 

-deg  F    R] 

?M  =  1000 

T     4. 

sat 

=  100°F 

Q 

=  9678.9  Btu/hr 
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TABLE  16 


Temperature  Distribution  Within  the  Fin  at  the 


10 


th 


Increment 


(Stainless  Steel  Condenser  with  25  Finite  Elements) 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

99.79 

15 

95.17 

2 

99.08 

16 

94.44 

3 

99.20 

17 

93.94 

4 

99.19 

18 

93.91 

5 

98.24 

19 

93.83 

6 

98.37 

20 

93.69 

7 

97.23 

21 

93.11 

8 

97.21 

9 

97.10 

10 

96.66 

11 

95.48 

12 

95.50 

13 

95.47 

14 

95.36 

a  =  50 

h     4-   = 

out 

deg 

1000  But/1 

ir-ft2- 

T   =  70°I 

GO 

-deg  F    RPM  =  IOC 

1 
)0 

T     = 

/sat 

100°F 

-  Q  =  8186. 

3  Btu/hr 
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TABLE    17 


Temperature  Distribution  Within  the  Fin  at  the 

10"*1"  Increment 
(Stainless  Steel  Condenser  with  40  Finite  Elements' 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

99.82 

15 

96.70 

29 

93.62 

2 

99.31 

16 

95.73 

30 

93.26 

3 

99.39 

17 

95.50 

31 

93.12 

4 

98.67 

18 

95.50 

5 

98.72 

19 

95.48 

6 

98.84 

20 

95.41 

7 

97.97 

21 

95.30 

8 

97.99 

22 

95.12 

9 

98.03 

23 

94.62 

10 

98.  09 

24 

94.45 

11 

97.25 

25 

93.93 

12 

97.24 

26 

93.91 

13 

97.20 

27 

93.85 

14 

97.08 

28 

93.75 

a  =  50 

h    4.   = 

out 

deg 

1000  Btu/1 

ir-ft2- 

T 

< 

deg  F     R] 

=  70°F 

•a 

?M  =  10  00 

T   .  * 
sat 

■  i  ■■— 

100°F 

Q 

=  8164.6  Btu/hr 
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TABLE  18 


Temperature  Distribution  Within  the  Fin  at  the 


10 


th 


Increment 


Stainless  Steel  Condenser  with  108  Finite  Elements) 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

99.84 

24 

97.21 

50 

94.46 

2 

99.50 

26 

97.13 

51 

94.39 

3 

99.54 

28 

96.64 

52 

94.64 

4 

99.09 

29 

95.83 

58 

94.29 

6 

99.19 

30 

95.52 

60 

93.83 

7 

98.65 

31 

95.43 

61 

93.76 

8 

98.67 

32 

96.32 

62 

93.87 

10 

98.81 

34 

96.28 

63 

93.84 

11 

98.19 

36 

96.15 

64 

93.79 

13 

98.23 

38 

95.82 

65 

93.73 

15 

98.36 

39 

95.31 

66 

93.64 

16 

97.70 

40 

95.03 

67 

93.55 

18 

97.72 

41 

94.95 

68 

93.31 

19 

97.73 

42 

95.45 

69 

93.31 

i 

21 

97.73 

44 

95.41 

70 

93.13 

22 

97.22 

48 

95.04 

71 

93.07 

a  =  50 

deg 

T 

CO 

=  70°F 

h  4.  = 
out 

1000  Btu/] 

ir-ft2- 

-deg  F    RP 

M  =  10( 

30 

T      4.    = 

100°F 

-  Q 

=  8141 

6  Btu/hr 
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Figure  20   Heat  Transfer  Rate  (Q)  of  Internally 

Finned  Condenser  at  a  Particular  Value 
of  h    ,  vs.  Saturation  Temperature  of 


out 


Working  Fluid 
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RPM  =  1000  2 

h    =  1000  Btu/hr-ft  -deg  F 

=  copper 

-  -  -  -   stainless  steel 
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CO 
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a 


=  10 


a  =  30° 
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Figure  21 


Heat  Transfer  Rate  (Q)  of  Internally  Finned 
Condenser  at  a  Particular  Value  of  Fin  Half 
Angle  vs.  Saturation  Temperature  of  Working  Fluid 
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Figure  22   Heat  Transfer  Rate  (Q)  of  Internally- 
Finned  Copper  Condenser  vs.  RPM 
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Figure  23   Heat  Transfer  Rate  (Q)  of  Internally  Finned 

Stainless  Steel  Condenser  vs.  RPM 
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Figure  24   Thermal  Resistance  of  Internally  Finned 

Copper  Condenser  vs.  h 
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Figure  25   Thermal  Resistance  of  Internally  Finned 

Stainless  Steel  Condenser  vs.  h 
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Figure    27      Comparison   of    Heat    Transfer   Rate    (Q-p,/Qc) 
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were  varied.   In  Figure  20,  the  influence  of  the  outside 

convective  heat  transfer  coefficient,  h   . ,  is  shown  both 

out 

for  copper  and  stainless  steel.   Figure  21  illustrates  how 
the  rate  of  heat  transfer  was  affected  when  the  fin  half 
angle  was  varied.   In  addition,  the  heat  transfer  rate  was 
found  not  to  vary  significantly  with  rotational  speed  when 
the  heat  pipe  operated  with  low  values  of  outside  convective 
heat  transfer  coefficient,  h    .   However,  for  high  values 
of  h    ,  the  effect  of  rotational  speed  cannot  be  ignored 
(see  Figures  22  and  23) . 

The  small  variance  of  the  heat  transfer  rate  with  rota- 
tional speed  is  due  to  the  dominancy  of  the  outside  thermal 

resistance,  R   ,  ,  such  that  changes  in  the  condensate  film 
out  3 

thickness,  5,  will  not  significantly  affect  the  total  thermal 

resistance  as  shown  in  Figures  24  and  25.   The  heat  transfer 

rate  will  have  a  large  variance  when  the  heat  pipe  is  operated 

with  a  high  value  of  h 

3  out 

From  the  results  of  the  parametric  study,  as  shown  in 
Figures  27  and  28,  the  heat  transfer  rate  continuously  in- 
creases as  the  fin  half  angle  decreases  (b/a  increases). 
Since  the  total  number  of  the  fins  in  the  condenser  is 
larger  with  a  small  half  angle  fin,  this  results  in  an  in- 
crease in  the  heat  transfer  rate.   However,  the  increase 
is  only  slight  when  the  fin  half  angle  less  than  11  degrees 
(b/a  greater  than  5) . 

The  study  shows  the  enhancement  in  heat  transfer  is 
larger  when  operating  with  a  high  value  of  hQut  or   when 
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R    is  small.   Therefore  the  design  operating  with  a 
high  value  of  h    ,  using  both  internal  and  external  fins 
is  recommended. 
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IV.   CONCLUSIONS 

1.  The  Finite  Element  Method  works  and  converges, 

2.  The  designer  has  a  choice  of: 

-  material 

-  number  of  fins  and  fin  angle 

-  RPM 

~   h     4- 

out 

3.  The  heat  transfer  rate  increases  continuously 

with  decreasing  fin  angle;  however,  for  half 
angle  less  than  11  degrees,  the  increase  is 
only  slight. 
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V.   RECOMMENDATIONS 

1.  Manufacture  internally  finned  heat  pipe  and  test. 

2.  Manufacture  internally  and  externally  finned  heat 
pipe  and  test. 
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APPENDIX 
COMPUTER    PROGRAM    LISTING 


**************  ***************************************^ 

* 

* 
TWO  DIMENSIONAL  HEAT  CONDUCTION  FINITE  ELEMENT  ANALYSIS  * 
OF  ROTATING  HEAT  PIPE  ,  USING  TRIANGULAR  ELEMENT  MODEL  * 
COMPILED  BY  MAJOR  IGNATIUS .S . PURNOMO  IN  JUNE  ,1973        * 

* 
* 
****** t  *******  l;***:K************************************x***** 
IMPLICIT  REAL*3(A-H,0-Z) 
DIMENSION  Z' 200  )  -  EPS ( 20U ) ,  HZ (200 ) , XC0FC5 ) , COF< 5 ) . R00TR(4  ) , ROOT I ( * 

%.K  Z^O ) . GE { 200 ) ■ DMDOT ( 200 ) , UF ( £00 ) , CF( 200 ) , RHOF ( 200 ) , DEL  C 200 ) . CU~ 

*00 )  •  AMTOT  ?  200  ) . R f 200  )  ■ QINC ( 200 ) * TB< 200 ) , TTC 200 ) , TBM ( 200 ) ,  TE  ( 200  )~ 

*IBC20C ) .NCC2O01 
DIMENSION  TSS< 3  '  . FANGLX3  )  ,HINF(3 ),MRPM(3) 

COMMON/  APOL/DOBF  ,  DOTH, KFINi  SO  ) ,  KFF  ( 50  ) ,  IFF  .  JTC ,  JLC ,  J  INT ,  KT 
COMMON/MAFO/A (  ^00 ,  50 ) ,  F  ( 200  .  1  ) ,  H  ( 200  ) *  TS<  200  > .  TSAT ,  CK ,  NEL ,  NSNP  .  M^ 

*N. ICOR(2O0.3) 
COMMON/PCRD'X  ( 20O  ) ,  Yf  20O  ) ,  EZEPO ,  BUIN ,  THICK  ,  ThLFA  ,  APS 
HDEN(Al,Bl,ZZ)-(-l  .OD0*HAl*ZZ**3-'30D0+Bl*ZZ**2/2  ODO  >  ) 

ELEMENT  CONNECT I U IT I ES 

READCS, 100)  NEL, NSNP, MB AH 

FORMAT (315) 

URITEC6.1S0)  NEL, NSNP. NB AN 

FORMAT (/2X. 'NO . OF. ELEMENTS- '. 15* 10X, 'NO. OF  SYSTEM  M.P-' . I5.10X.  - 

*0  OF  BANDED*' , IS) 
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READ (S. 200)  i  IEL. CICORf IELj I), 1-1,3). IEL-l.NEL) 

F0RMmT(4I5) 

URITEC6.2S0) 

FORMAT*. -'2X-  'ELEMENT'  ,  1CX.  'NP1',14X,  'NP2'  .  15X,  'NP3'  ) 

URITE'b/251  KIEL,  CICORCIEL.  I),  1-1.3).  IEL-l.NEL) 

FORMATCIS, 12X, 15, 12X , 15. 12X.I5) 

THE  CONDENSER  GEOMETRY 

CL  IS  THE  CONDENSER  LENGTH 

CfflGL    IS    THE  CONE  HALF  ANGLE 

RBASE  IS  THE  CONDENSER  RADIUS  AT  THE  BASE 

THICK  IS  THE  UALL  THICKNESS  OF  THE  CONDENSER 

READ < 5 • 300  )  CL , CANGL . RBASE , R2 , THICK , BFIN 

FORMAT' 6G10  5  ) 

CL  =  CL-'12  ODO 

R2  =  R2^'12  ODO' 

RBA5E-RBASE/12  ODO 

B"IN=BFIN/12  ODO 

'  NDIU  IS  THE  NUMBER  OF  INCREMENT 
NEST  IS  THE  ELEMENT  NUMBER  AT  THE  END  OF  THE  TROUGH 
NEFB  IS  THE  ELEMENT  NUMBER  AT  THE  BASE  OF    THE  FIN 
NBOTI  IS  THE  FIRST  ELEMENT  AT  THE  BOTTOM  SIDE 
NBOTF  IS  THE  LAST  ELEMENT  AT  THE  BOTTOM  SIDE 

READ( 5 . 400  )  NDIU -  NEST , NEFB , NBOTI , NBOTF 

F0RMATC5I5) 

DIU  =  DFLOAT(NDIU  ) 

PI-3  1415926S358979D0 

PHI =2  0D0*CANGL*PI/3c0  ODO 

SPHI-DSINCPHI ) 

CPHI»DCOS(PHI ) 

TPHI-DTANCPHI ) 
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DELX-CL'Dlu 
C£ASE*2. 0D0*PI*RBASE 
REXIT*RBASE+CL*TPHI 
CEX I T -2 . 0D0*P I  *REX  1 T 
THICK-THICK/12.0DO 

URITE<6,450)CL,CANGL.R2.THICK.RBASE,BVIN 
50   FORMATS     \/V5X,  'CONDENSERS    LENGTH- '  ,£12. 5, 30X,  'HALF   CONE   ANGLE- 

*,E12   5,/.SX, 'R2-' ,E12   S,45X, 'THICK- ' , E12   5, / , SX, 'R-BASE- ' ,E12 . 5. - 

*X, 'B-FIN«',E12   S) 

DATA  FOR  RUNNING 

READ' 5- 500)  (MRPMCI  ).  1-1,3) 
'50  FORMAT'  315)  ■ 

READ (5, 60O)  (TS5CI  ), 1-1.3) 
i30  FORMAT  (3G10  5) 

REaO<5.7O0)  TINF 
00  FORMATCG10  5) 

READ (5,300)  (HINFCI  ). 1-1,3) 
30  FORMAT (3G10. 5) 

THE  CONVERGENCE  CRITERIAN 

READ(5,900)  CRIT 
30  FORMAT (GIO  9) 

INTERNALLY  FIN  GEOMETRY 

READ(5,1000)  (FANGLCI  ), I»l,3)  . 
30  FORMAT C3G1 0.5) 

READC5.1100)  IFF 
50  FORMAT(IS) 
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READ'S, 1200)  (KFIH< I  )  JCFFCI ;,I-1  -IFF) 
20  F0RnAT(16IS) 

READ  •;  5 ,  1300  iDOBF  ,  DOTH  .  JTC  ,  JLC ,  JINT ,  K7 
30  F0PriHTiaGi0.5.4I5) 
,  NHB-NEFB/2 

'nTN-NBOTI+CNBOTF-NBOTI >/2 
NEF-NBOTF+1 
URITEC6, 1350)  ICDRCNBOTI . 2  ) , IGOR (KERB, 1 ) . ICORCNTtt, 2 ) . ICORCNEST. 1~ 

XICGRCNBOTF-l  ) 
I'd   F0RflAT<///5X,  'TIB=y  , I5,10X,  'TT*M5,/.5X,  'TBW-',  15-  10X,  'TE-M5,a 

*6X. 'TB«' . 15) 
DO  1400  KIA«1,3 

ALFA=FRNGL(KIA  >*2  0D0*PI/360  0DO 
SALFA ■ DS INC  ALFA ) 
CALFA=DCOSCALFA) 
TALFa-DTAHihLFA.) 
EZER0=>2  ODO*B'JIN*TALFA 
EPSO=EZERO 

NFIN-CBHSE-'  ( EZEPO+EPSO  ) 
EPSEX-fCEXIT-(NFIN*EZERO) )/NFIN 
BETA»  * EPSEX-EPSO  J/DIU 
ZZERO-BUIN/CALFH 
ZA=0  OD0 

BOUNDARY  CONDITIONS  AND  TEMPERATURE'S  ESTIMATE 
ALONG  THE  FIN  BOUNDARY 

DO  1450  NTINF-NBOTI,NBOTF 
3  TSCNTINF)«TINF 

DO  1500  NNT-NSF.  MEL- 
TS (NNT)-O  0D9 
<9  H(NMT)»0  0D0 
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DO  1600  IGT-1,HEST 

IE-IC0R(IGT,2) 

READ'S^ 1700)  T(I£) 

IG-IC0R(NEST, 1  ) 

PEhO'S,  1700)  T(IG) 

FOPNATcGlO .5) 

DO  mOO  KIB»1,3 

NRPN*riP.PM<KIB  ) 

0MEGA»NRPnt2    0DO*PI*60   ODO 

DO    1400   KIC-1,3 

DO    1200   KL-NBOTI,NBOTF 

H(KL)»HIflFCKIC  ) 

HIFN-HINFCKIC > 

DO    1400   KID«1.3 

TSHT-TSSCKID > 

DO  1900  NSAT=1,NEST 

TSCNSAT  )*TSAT 

TSOLID= (T3A7+TINF  )/2 . ODO 

DEL(l>-0  0O006752D0 

QT-J  ODO 

QTOT=0  ODO 

DHTOT-O  ODO 

NK-NDIU+1 

DO  2000  NI-1,NK 

R(NI  )»R2+NI*DELX*SPHI 

EPSCNI  )«EPSO+NIXBETA 

APS°EPS(NI j 

NODAL  POINTS  COORDINATE 

CALL  COORD 

2(1 )«ZA 

DO  2100  IZEL-l.NEFB 

NA-ICORCIZEL. 1 ) 
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NB-ICORaZEL.2) 
XE-X<NA  >-X(NB ) 
VE«V'.NA)-Y<NB  ) 
ELZ«DSGPTO<E**2+YE**2 ' 
j  ZaZEL+1  >=Z<  IZED+ELZ 
XZB-XLlCORCNHB.l  )  )-X(  ICOR(  1 .2  )  ) 
VZB«V«'IC0R<MH&,  1  i  )-Y<  IG0RC1.2) ) 
ZE=DSQPT< XZB**2+YZB**2  ) 
ZOZ2EP0 
IM-1 

PARABOLIC  TEMPERATURE  DISTRIBUTION  ALONG  THE  FIN 
BOUNDARY. USING  LAGRANGE  INTERPOLATION 

TPl-Tr I COR  CI, 2 >  ) 

TP2*T(  ICORCNHB. 1  I  ) 

TP3  =  TriC0R(NEFB, 1  )  ) 

AP1-TP1/(ZB*ZC  ) 

fiP2  =  TP2/C2B^»:ZB-2C)  ) 

AP3=TP3/(2CX(2C-2B)  ) 

BPl  =  -iZB  +  ZCUAPl 

BP2=-2C*hP2 

BP3»-2B**P3 

A1*AP1+AP2+AP3 

B1-BP1+BP2+BP3 

TC  =  0  0D0 

DO  2200  NY-1 .NEST 

TC-TC+T(IC0R(NY,2)  ) 

AY»DFLOATCNY+l  ) 

TF-  CTC+T( ICCRCNY.  1  )  )+AY*TS(NY  )  UiZ   0D0*AY  ) 

SOLID-FLUID   PROPERTIES 

HFG-1097.2D0-0    601875DO*TS( 1 ) 
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RH0F'NI>*62.774D0-0  0O2S5638DO*TF-O . 000053572D0*TF**3 

CF'.III  )=O.3034D0+0.00073S927D0*TF-0.00000147321D0#TF**2 

UFCNI  )-O.O01397D0-0  . 00001 4b69P0*TF+0 .0000000631253D0*TF**2-0,000~ 

*00Q0097G569D0*TF**3 

UFCNI  )-3600*UF(HI ) 

CUCNI )=231 . 7772D0-0 . O2222D0*TSOLID 

CK-CUCMI  ) 

CONST -RHCFCNI  )**2*0f1EGA**2*HFG*CPHI*CALFAXR(NI  ) 

AVERAGE  ELEMENT  CONUECTIUE  COEFFICIENT  ALONG  THE 
FIN  BOUNDARY 

ZSTAR»ZZERO-DEL(NI  )/CALFA 

AZZ-DELCNI  )/SALFA 

ZZ=ZST*R 

AZS»DABSC4*CF(NI  )*UF(NI  )*HDEN( Al , Bi>ZZ )/CONST )**0  25D0 

HAC=0  ODO 

DO  2300  IEL»1,NEFB 

AZ-Zf  IED 

BZ-ZCIEL+1  ) 

IFCZSTAR  LE.BZ)  GO  TO  2400 

GO  TO  2500 
!-)0  IF  (MAC  NE  0  ODO)  GO  TO  2600 

BZ-ZSTAR 
|50  IFCIEL  NE  1)  GO  TO  2700 

IAK»(BZ-AZVS.0D0 
ZZ-AK 
GO  TO  2SO0 
!"30  A<-(BZ-AZ)/4.OD0 

ZZ-AZ 
:iJ0  ZEL=4*AK 

DO  2900  NH-l.S 

HZ(NH)»DABSCCF(NI  )**3*C0NST/(4*UF(NI )*HDEN( Al , Bl , ZZ ) ) )X*0 . 25D0 
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i   ZZ-ZZ+AK 
CONH«AiaCHZ(l  >+4#HZ(2)+S*HZ<3)+4*HZ<4)+HZ(Sn/<3*ZEL) 
IFCZSTAR    EQ.BZ)    GO   TO   3000 
H(IEL)*CONH 

GO   TO   23M 
i    AZ»ZSTAR 
HAZ«CONH*<  AZ-Z'.  IZL)  ) 
DELA*AZ3 

BZ»Z(IEL+1 ) 

DELB=  C BZ-ZSTAR  )*AZZ/(ZZERO~ZSTAR ) 

DELZ-  (.  DELA+DELB  >/E  ODO 

HAC=« ( BZ-AZ  )*CF ( NI )/DELZ 

H( IED-CHAZ+HAC  )/(EZ-Z(IEL) ) 

GO    TO    5300 

AZ-ZCIEL1 

DELA-DELB 

HAZ=0.0DO 

GO   TO    3100 

CONTINUE 

NETI-NEFB+1 

DO    3E00    IEL-NETI.NEST 

H(IEL)-CF(NI J/DELCNI ) 

ENTRY  INTO  THE  FINITE  ELEHENT  SOLUTION 

CALL  FORMAF 

CALL  BANDEC ( NSNP . NBAN . 1 ) 

THE  TEHPERATURE  DISTRIBUTION 


DO  3300  NT- 1 -NSNP 
l!$0  T(NT)«F(NT,  1  ) 

TIBCNI  )-T(ICCR(NB0TI.2)) 
TTCNI  )«T(ICOR(NEFB< 1  )  ) 


99 


TBNCNI  )-T(IC0P(NTM,2J  ) 
TE(NI  )=Ti  ICOR'  riEST,  1  )  ) 
TE<MI  '-T<ICOR<NBOTF,l ) ) 
TTS-0    OD0 
DO   3400   NS-l.NSNP 
bo   TTS-TTS+T(NS) 
PN-DFLOATCNS  > 
TS0UD-TT3/PN 

Q  AT  THE  BOTTOM  SIDE 

QBI'O  OD^i 

DO    3500    IEEL=NBOTI,HBOTF 
NKA-ICORdBEL,  1; 
NKB«IC0RCIBEL,2) 
XB-XCNKA  >-X(NKB ) 
YB-VCNKA J-V(NKB ) 
ELB=DSQRT(XBX*2+Y3**2  ) 
50   QBI«GEH-CT(NKA)+T(NKB)-2*T3(IBEL)  )*ELB*H<  IBED/2  .  0D0 
QB(NI  )«QBI*DELX 

ITERATION   UNTIL   CONVERGENCE   CRITERIA    IS   MET 

IFCiri   EQ    1  )    GO   TO   3600 
QJ-QBI 
GO    TO    3700 
3130   QI-QBI 

in»2 

GO  TO  2 
T30  AQ-DA3SCQJ-QI )/QJ 

IFCAQ  LE  CRIT)  GO  TO  3800 

QI-GJ 

GO  TO  2 
m  DHDOTCNI  )«2*QBlXDEtX/HFG 
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0 

00SDO 

AMD 

ROCTR ( 1  ) 

GT 

0 

0D0) 

GO 

TO 

39O0 

o 

0O5D0 

AMD 

R00TRC2) 

GT 

0 

0D0) 

GO 

TO 

4000 

0 

005DO 

•,ND 

RCCTRC3) 

GT 

0 

OD0) 

GO 

TO 

4100 

0 

005D0 

AND 

R00TRC4) 

GT 

0 

OD0) 

GO 

TO 

4800 

DfTTOT-DMTOT+DnDOU  MI  > 

Cl«RHOF(MI )**2*0MEGA**2*R(NI )/(3*UF(NI ) ) 

XCOFC1  >«-DMTOT 

XC0FC2  »-0.  OD0 

XCOF<3)-0   0D0 

XCOFM  »»C1*EPS(MI  ) 

XC0F15)*C1*TALFA 

MM 

CALL  DPOLPT ( KCOF , COF , M , ROOTR , ROOTI ,  IER  ) 

IF < ROOTR (1  )  .LT 

IFCROOTRra)  . LT 

IF (ROOTR (3  )    LT 

IF (ROOTR (4)  LT 

URITEfb  4300) 
>0  FORMAT* //10X,  'CPASH, CRASH, CRASH'  ) 

URITEC6.4350)  ( ROOTR ( I ) , 1-1 , 4 ) 
J0  F0RriAT(//5X.4(E12.7,3X)) 

GO  TO  6100 

THE   CONDENSATE   THICKNESS 

>0   DELCNI+1  )- ROOTR U  ) 

GO    TO    4400 
10    DELCNI  +  1  J-R00TRC2) 

GO    TO    4400 
>0   DEL C  NI  + 1  J -ROOTR ( 3 ) 

GO   TO    4400 
)0   DELCNI  +  1  ;=R00TR(4) 
10    QEL  =  0   0D0 

IFCNI    NE    10)    GO    TO    4500 

URITECS.4600) 
>0   FORHATC  'l',//5X,  'NP' ,6X,  'X'  , 12X. 'V . 12X,  'T'  ) 

DO    4700    NP-l.NSNP 
.0    URITEC6.4300)    NP.XCNP). V(NP) ,T(NP) 


101 


FORflAT  ( /2X  ,13,  3X .  3 (F 10  . 6 ,  3/! )  ) 

URIT£«'6,4900) 

F0RMAT<//2X,  'EL',8X,  'H',UX,  ' EL-LENGTH ',  15X ,  'Q-EL'  ) 

DO  5000  KKL-l.MBOTF 

hkx»icor(kkl.  1 ) 

nky-icor(ki:l.2) 

XP-X<r{KX)-X(NKY) 
VP-Y(NKX)-YCNKV) 

EXV«DSQRTf  XP**2+VPt *2 ) 

QEP-DrtBS r  ( T ( NKX  ^+T ( NKY  )-2XTS ( KKL  )  )*EXY*H<  KKL  >/2  .  OD0  ) 

QEP«QEP*DELX 
fed  URITE(6,510O>  KKL  .  H<  KKL  • ,  EXY  .  OEP 
DO  F0RPIAT(/3X.  12. 3X, E12. 5, 3X.E12  5.10X.E12.5) 

URITEf.6-5200)  CRIT 
DO  FORMAT (/2X.  'CONUERGENCE  CRITERIAN- '  ,E15  .  8 ) 

Q  FROM  THE  TOP  SIDE 

50   DO   530O    IQEL-l.NEST 

KA-ICORCIGEL. 1 » 

KB-IC0RCIGEL,2) 

XQEL-XCKA  >-X(KB) 

YOEL-YCKB  )-V(KA) 

ELM-DSQRT(XQEL#*2+YQEL**2) 

QEL-QEL+(2*TS(IQEL )-T(KA  )-T(KB ) )*ELM*H< IQEL )/2 . ODO 
)0  CONTINUE 

QINCCNI  )-QEL*DELX 

AHTOTCNI  )-DffTOT 

GET-QEL*DELX*NFIN*2 

QT-QT+QET 

QA-QBI*DELX*NFIN*2 

QTOT-QTCT+GA 

CONTINUE 

URITE £ G . 5400  )  HFG . NF IN . H ( NBOTF  ) , T5AT , HRPH , QTOT , QT , FANGL ( K I A ) 
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> 

;4'0  FORMAT*.  '  '  »//,5X,  'HFG-'.E12  5,  /,  SX,  'NUMBER  OF  FIN»  '  .  15,  /,  SX,  'H-0~ 

*•  '  .E12  S,/,5X,  'TSAT»'  ,E12.5,/,5X.  'RPM-'  , IS,/,5X,  'G-BOT-'  , E12.5./~ 
>  i 

XX, 'Q-TOP-' ,E12 5,/,5X. 'HALF-ANGLE-  '  ,  F3  .  3  ) 
WRITE (6,5500) 
5E'0  FORMAT (  '0',6X,  M\4X,  'FILM  THICKNESS ',  6X,  'Q-INCREM'  ,  6X,  'MASS-TOT^ 

17X, 'TIP' ,3X. 'TT' . 10X, 'TE'.2X. 'TB' ) 
DO  5600  NR-1 -NDIU 
"6'0  UPITE  <  6  •  555  )  MP ,  DEL  <  NR  ) . OB ( NR  ) ,  AMTOT  CNR ) ,  TIBCNR ) ,  TT  CNR  ) ,  TECNR ) ,  T* 

XNP  ) 
;=;   FORMAT  ('  ',4X,I4,4X,F12  10,  4X,  F10  .  4,  6X,  F9  5,  6X,  F5  . 1 ,  6X,  F5  .  1 ,6X,  /> 

1 . 1,6X,F5  1 ) 
UR I TEf 6.5SC0  5 
:S>0  FORMAT (  '0'  .SX,  'J'  ,6X,  'K-UALL'  ,  4X,  'K-FIUV  ,  3X,  'DENSITY'  ,  4X,  'UISC-~ 

XILM' .6X.  'EP5IL0N' ,SX, 'RADIUS' ,5X, 'TBM' ,5X,  'Q-BOT' ) 
DO  5900  NG-1,NDIU,2 
:S)0  URITEC6.777)  NG.  CU(NG  ) ,  CFCNG  )  ,RHOF(NG  ) ,  UF(NG  ).  EPSCNG  ) ,  RCNG  ) ,  TBflC  ~ 

X)  -QINCCNG) 
"'   FORMRT  ('  ' ,4X. I4.4X.F7.3.4X,F6 4,4X,F6  3, 4X, F9 . 7, 4X. F9 . 7, 4X, F7  - 

15X,F5  1.1X.1P1D12.3) 
bo  CONTINUE 

;:.J0  stop 

END 

SUBROUTINE  COORD 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/PCRD/X ( 200  )  • V C 200  )  - EZERO . BUIN . THICK , TALFA . APS 

COMMCN/APOL/DOBF , DOTH , KFIM C 50  ) . KFF ( 50 > , IFF . JTC , JLC . JINT . KT 
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DELH-BUIH/DOBF 

X(1)=0  OD0 

Yi.  1  )»THICK+BUIN 

N-i 

DO  321  1*1. IFF 

ZCA-KFIHCI  ) 

ICB-KFF"  I  ) 

CEA=DFLOAT< ICB-ICA ) 

AN-0  0D0 

DO  322  II-ICA, ICE 

X(II )«=X< 1 )+N*AN*DELH*TALFA/CBA 

V(II )=V(1  )~N*D£LH 

AN-AN+1 . 0DO 

N»M+1 

AN-0  0D0 

ICD-ICB-ICA+1 

DO  323  J=JTC, JLC.JINT  • 

X(J )-X( 1  ) 

V(J)-<1  OD0-AN/DOTH)*THICK 

DO  324  JJ-1. ICD 

X  ( J+JJ  )  -X  (  J  )+ J JXEZERO-'  ( 2*  ( CBA+1 .  OD0  )  ) 

Y(J+JJ)«Y(J) 

DO  335  K«1,KT 

X C  J+JJ+K J -X ( J+JJ  )+KXAPS/ ( 2 . 0D0XKT ) 

V(J+JJ+K)-Y(J) 

AN«AN+1  0DO 

RETURN 

END 

SUBROUTINE  FORHAF 

IMPLICIT  REAL*8(A-H,0-Z) 

DIHEMSION  B(3),C(3).EA(3.3) 

connoN/PCRD/x < 200 ) . y c 200 » . ezero , buin . thick , talfa , aps 

COrinON/MAFO/A ( 200 . 50  ) . F ( 200 , 1 ) , H( 200 ) . TS ( 200 ) , TSAT . CK . NEL . NSNP , N~ 
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*H, I COR (500. 3) 
DO  101  N-1,NSNP 
F(M, 1  i-0  ODO 
DO  102  MA -1, MEAN 
A(N.MA)«O.0D0 
CONTINUE 

DO  103  IEL-1,NEL 
IA-ICOR'  IEL,1  i 
IB-ICORv IEL,2) 
IC-IC0RCIEL.3) 
B(l  )-Vv IB  l-V(IO) 
E(2)=V(IC1-VC IA) 
Bc3  i-vartJ-VdB) 
C(l  )«X(IC  >-)<(  IB  ) 
C<2)-X(. IA J-XCIC  ) 
C(3)-X(IB)-XCIA) 
EL-DS0PTCCC3  )**2+B< 3 )**2) 
AS-DABSC (B(l  )*C(2)-B(2)*C(1 ) 1/2.0D0) 
HC-HC  IED/CK 
DO  104  J-1,3 
JJ-ICORCIEL,J) 
DO  105  K-1.3 
KK«ICOR(IEL,K) 

EA <  J  ,  K ) - ( B ( J ) *B ( K)  +C ( J  )  *C ( K ) ) / ( 4XAS > 
IFCHC  EG  0  0D0)  GO  TO  106 
HEL«HC*EL/6  OD0 
IF(J  EQ  3)  GO  TO  10S 
IFCK.EQ  3)  GO  TO  106 
IFCJ.EQ  K)  GO  TO  107 
EA(J.K)»EA(J.K)+HEL 
GO  TO  106 

EA(J,K)»EA(J,IO+2*HEL 
IF(KK.LT.JJ)  GO  TO  105 
NU-KK-JJ+1 
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A  ( J  J ,  NU  '  -  A  ( J  J ,  NU  )  +EA  ( J  •  K  ) 

CONTINUE 

CONTINUE 

FE»HC*TS(IEL.UEL/2  ODO 

Fa  A.  1  >=F<IA.  1  )+FE 

F(ID,  1  >*FCIE,  1  i+FE 

CONTINUE 

RETURN 

END 

SUBROUTINE  BANDEC ( NEQ . NAXB - NUEC ) 

IMPLICIT  REAL*3CA-H,0-Z> 

COMMON-  PCRD'  X(  £00  ) ,  V(20tf  ) .  EZERO,  BUIN,  THICK  .  TALFA.  APS 

COnnON/NrtFO/  A  <  200 , 50  ) . F ( 200 ,  1  ) ,  H<  200  ) ,  TS  < 200  ) ,  TSAT ,  CK ,  NEL ,  NSNP ,  N^ 

*N, ICORC2^0, 3  > 
LQOP-NEQ-1 

DO  201  I -1, LOOP 

MB-I+1 

riB«MIN0<  I+MAXB-1  >  NEQ  ; 

DO  201  J -MB. MB 

L-J+E-MB 

D-A<I,L)/A(I, 1 ) 

DO  202  Hn-i.NUEC 

f  c  j .  nn )  ■  f  (j ,  nn  >-d*f  c i . mm  ) 
mm-mimocmaxb-l-h.neq-j+i ) 

DO  201  K-l.MM 

NN-L+K-1 

A<  J,K)-A( J.K  )-D*A<I,NN) 

DO  203  1-1 ,  NUEC 

FCNEQ. I  )-F(hEQ. I  )/A(NEQ.l ) 

DO  204  I =2. NEQ 

J-NEQ-I+1 

K'MINOCNEG-J+l-MAXB) 

DO  204  mm«i,nuec 
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DO  205   L  =  2,K! 

riE-j+L-i 

F<  J  ■  MM  i-F ( J>  (in  '-h(  J .  L  >*F(  MB ,  MN  ) 

fj.fin  j-Fcj.rinj/'ACj,  1 ) 

RETURN 
EHD 
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